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ABSTRACT 

We present detailed radiative transfer spectral synthesis models for the Iron 
Low Ionization Broad Absorption Line (FeLoBAL) active galactic nuclei (AGN) 
FIRST J121442.3+280329 and ISO J005645. 1-273816. Detailed NLTE spectral 
synthesis with a spherically symmetric outflow reproduces the observed spectra 
very well across a large wavelength range. While exact spherical symmetry is 
probably not required, our model fits are of high quality and thus very large 
covering fractions are strongly implied by our results. We constrain the kinetic 
energy and mass in the ejecta and discuss their implications on the accretion 
rate. Our results support the idea that FeLoBALs may be an evolutionary stage 
in the development of more "ordinary" QSOs. 

Subject headings: AGN: FIRST J121442.3+280329, ISO J005645. 1-273816 

1. Introduction 

Spectroscopic observations of quasars show that about 10-20% have broad absorption 
troughs in their rest-frame UV spectra (see Trump et al. 2006, for example). These ab- 
sorption lines are almost exclusively blueshifted from the rest wavelength of the associated 
atomic transition, indicating the presence of an outflowing wind in our line of sight to the 
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nucleus. The line-of-sight velocities range from zero to up to tens of thousands of kilometers 
per second (e.g., Narayanan et al. 2004). 

While understanding these outflows is of fundamental interest for understanding the 
quasar central engine, it is also potentially important for understanding the role of quasars 
in the Universe. The observation that the black hole mass is correlated with the velocity 
dispersion of stars in the host galaxy bulge (e.g., Magorrian et al. 1998; Ferrarese & Merritt 
2000; Gebhardt et al. 2000) indicates a co-evolution of the galaxy and its central black hole. 
The close co-evolution implies there must be feedback between the quasar and the host 
galaxy, even though the sphere of gravitational influence of the black hole is much smaller 
than the galaxy. Energy arguments, however, show that is quite feasible that the black hole 
can influence the galaxy; as discussed by Begelman (2003), the accretion energy of the black 
hole easily exceeds the binding energy of the host galaxy's bulge. 

The nature of the feedback mechanism that carries the accretion energy to the galaxy 
is not known. Since AGNs are observed to release matter and kinetic energy into their 
environment via outflows, it is plausible that these outflows contribute to the feedback in 
an important way. One of the difficulties in using quasar outflows in this context is that 
they are sufficiently poorly understood that there are significant uncertainties in such basic 
properties as the total mass outflow rate and the total kinetic energy. 

What is the kinetic luminosity of the broad absorption line quasar winds? That turns 
out to be very difficult to constrain. While the presence of the blueshifted absorption lines 
unequivocally indicates the presence of high- velocity outflowing gas, the other fundamentally 
important properties of the gas, including the density, column density, and covering fraction 
are very difficult to constrain. 

The density is difficult to constrain because the absorption lines are predominately 
resonance transitions, and their strengths are not very sensitive to density. Without knowing 
the density, the distance of the gas from the central engine cannot be constrained; the same 
ionization state can be attained by dense gas close to the central engine, or rare gas far 
from the central engine. Density estimates are possible when absorption lines are seen from 
non-resonance transitions, but even then, they can differ enormously. For example, de Kool 
et al. (2001) analyzed metastable Fe II absorption lines in FBQS 0840+3633, and inferred 
a electron density < 1000 — 3000 cur 3 and a distance from the central engine of several 
hundred pc. In contrast, Eracleous et al. (2003) analyze the metastable Fe II absorption 
in Arp 102B with photoionization models and infer a density of at least 10 11 cm -3 and a 
distance of less than 7 x 10 16 cm. 

The global covering fraction is also difficult to constrain directly from the quasar spec- 
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trum; we know the gas, at least, partially covers our line of sight, but we have little infor- 
mation about other lines of sight. Covering fraction constraints are generally made based 
on population statistics. In a seminal paper, Weymann et al. (1991) showed that for most 
BALQSOs, the emission line properties are remarkably similar to non-BALQSOs. Thus, the 
fact that 10-20% of quasar spectra contain broad absorption lines is interpreted as evidence 
that there is a wind that covers 10-20% of sight lines to all similar quasars, and whether or 
not we see absorption lines depends on our orientation. Alternatively, some BAL quasars 
have notably different line emission than the average quasar; examples are the low-ionization 
BALQSOs studied by e.g., Boroson & Meyers (1992). These objects may instead represent 
an evolutionary stage of quasars, as the quasar emerges from the cloud of gas and dust in 
which it formed (Becker et al. 1997). 

While it seems that the column density should be easy to constrain, more recent work 
has shown that it can be very difficult to measure. It was originally thought that non- 
black absorption troughs indicated a relatively low column density for the absorbing gas 
(equivalent hydrogen column densities of 10 19 ~ 20 cm~ 2 , e.g., Hamann 1998). But it has now 
been found that the non-black troughs indicate velocity-dependent partial covering, where 
the absorption covers part of the emission region, and the uncovered part fills in the trough 
partially (e.g., Arav et al. 1999). Thus, the column density appears to be high, but it is very 
difficult to constrain directly from the data except in a few very specialized cases (see for 
example Gabel et al. 2006; Arav et al. 2005). 

How can we make progress on this problem? It is becoming clear that because of the 
difficulties described above, the traditional techniques for analysis of troughs (e.g., curve of 
growth) and modeling (e.g., photoionization modeling to produce absorption line ratios and 
equivalent widths) are limited. An approach that may be profitable is to construct a physical 
model for the outflow, and constrain the parameters of the model using the data. 

Our first foray into constructing physical models for quasar winds was performed by 
Branch et al. (2002). In that paper, the FeLoBAL 1 FIRST J121442+280329 was modeled 
using SYNOW, a parameterized, spherically-symmetric, resonant-scattering, synthetic spec- 
trum code more typically used to model supernovae (Fisher 2000). The difference between 
this treatment and a more typical one applied to the same data by de Kool et al. (2002) is that 
SYNOW assumes that emission and absorption are produced in the same outflowing gas. In 
contrast, the approach taken by de Kool et al. (2002) assumes that absorption is imprinted 
upon a typical continuum+emission line quasar spectrum; that is, the absorbing gas is sep- 
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arated from the emission-line region. In fact, based on the analysis of the Fe II metastable 
absorption lines, they find that the absorber is 1-30 parsecs from the central engine, much 
farther than the quasar broad emission-line region. Note that FIRST J121442+280329 is 
not the only object that can be modeled using SYNOW; Casebeer et al. (2004) present a 
SYNOW model of another FeLoBAL, ISO J005645. 1-273816. 

The SYNOW model is attractive because it is simple; only one component is needed to 
model both the emission and absorption lines. However, this model is limited. The primary 
purpose of the SYNOW program is to identify lines in complicated supernova spectra. Thus, 
individual ions can be added to a SYNOW run at will in order to see if features from emission 
and absorption from those ions is present. It does not solve the physics of the gas, so physical 
parameters beyond the existence of a particular species and its velocity extent cannot be 
extracted from the results. 

We test the ideas of Branch et al. (2002) and Casebeer et al. (2004) by using the 
generalized stellar atmosphere code PHOENIX to model the spectra of the two FeLoBALs that 
were successfully modeled using SYNOW, and including spectra that extend to rest-frame 
optical wavelengths for FIRST J121442+280329. PHOENIX is a much different code than 
SYNOW in that it contains all the relevant physics to determine the spectrum of outflowing 
gas. It solves the fully relativistic NLTE radiative transfer problem including the effects 
of both lines and continua in moving flows. For a discussion of the use of both SYNOW and 
PHOENIX in the context of modeling supernovae spectra, see Branch, Baron, & Jeffery (2003). 
We find that PHOENIX is able to model the spectra from these objects surprisingly well, and 
we are able to derive several important physical parameters from the model. 

In §2 we describe the PHOENIX model in detail. In §3 we describe our determination of 
the best-fitting model. In §4 we describe the results of our model fitting. In §5 we discuss 
the physical implications of the model, how it relates to other BAL spectra and where it fits 
in the BAL picture. An appendix includes a flowchart of a PHOENIX calculation. 

2. Models 

Photoionization codes have been essential in understanding emission and absorption 
features in the spectra of active galaxies. Cloudy (Ferland et al. 1998; Ferland 2003) is widely 
used. PHOENIX solves the radiative transfer equations exactly, and tries to have the most 
accurate radiative data possible. As a radiative transfer code, PHOENIX correctly produces 
line profiles due to relativistic differential expansion. It is this feature that we are making 
use of here. 
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In this paper we model the situation envisioned by Branch et al. (2002): the emission 
and absorption occur in a fairly optically thick expanding shell illuminated from the inside by 
the continuum; as discussed in §5.4, this situation may be a consequence of quasar evolution, 
occurring when the quasar ejects a shroud of dust and gas (e.g., Voit et al. 1993). 

PHOENIX is a spectral synthesis code; the direct output is a model spectrum. The only 
way to obtain fluxes or equivalent widths of lines in a PHOENIX model is to measure them 
directly from the synthetic spectrum in the same way that they are measured from the 
observed spectrum. Measuring emission and absorption lines from complex quasar spectra 
is well known to be rather uncertain, as a consequence of blending and uncertain placement 
of the continuum. So in PHOENIX, this step is bypassed, and the synthetic spectrum is 
compared directly with the observed spectrum. Second, the input parameters are somewhat 
different. In PHOENIX, density is given as a function of the radius in concordance with the 
assumed velocity profile as a function of radius. An analogy to the photoionizing flux is a 
little difficult to construct. As noted in the next section, two of the important parameters 
are the reference radius Rq, the radius at which the continuum optical depth at 5000A is 
unity, and the model temperature T mo d e i, defined in terms of the total bolometric luminosity 
in the observer's frame, L, and the reference radius. Thus, L or T are somewhat analogous 
to the photoionizing flux, because for a fixed reference radius, they give the intensity of 
the continuum at the reference radius. Finally, the column density can be evaluated for 
particular values of the optical depth. 



2.1. The Model Parameters 

Our models are spherically symmetric, with homologous expansion (v oc r). Homologous 
expansion is analogous to the Hubble expansion. The model atmospheres are characterized 
by the following parameters (see Baron et al. 2004, for details): (i) the reference radius 
i?o, the radius at which the continuum optical depth in extinction (r st( j) at 5000A is unity; 
(ii) the model temperature T model , defined by the luminosity, L and the reference radius, R , 
[r m odei = (-^/(47ri?o ")) 1 ^ 4 ]) where a is Stefan's constant; (iii) the density structure parameter 
v e , {p(v) oc e~ v / Ve ^]; (iv) the expansion velocity, v$, at the reference radius; (v) the pressure, 
P ou t, at the outer edge of the atmosphere; (vi) the LTE-line threshold ratio, equal to 5 x 10~ 6 ; 
(vii) the albedo for line scattering (metal lines only, here set to 0.95); (viii) the statistical 
velocity ( = 50 km s -1 , treated as depth-independent isotropic microturbulence, and (ix) 
the elemental abundances, assumed to be solar as given by Grevesse & Noels (1993). 

We emphasize that for extended model atmospheres one should not assign, a priori, 
a physical interpretation to the parameter combination of T modc i and Rq. While T mode i 
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has a well defined physical meaning for plane-parallel stellar atmospheres, its definition 
for extended atmospheres is connected to the particular definition of the radius Rq (see 
Baschek et al. 1991). In addition, the reference radius Rq in our models is defined using an 
extinction optical depth scale at A = 5000A and is not directly comparable to observationally 
derived radii. Therefore the model temperature is not well defined for extended atmospheres 
and must be regarded only as a convenient numerical parameter. We chose the maximum 
extinction optical depth so that model would be just optically thick to continuum scattering 
in order to replicate the optical spectrum observed in these objects (see § 4.1). 

2.2. Comparison of our model with de Kool et al. (2002) 

2.2.1. The de Kool et al. approach 

Because we want to compare and contrast our model with de Kool et al. (2002) we 
must first briefly describe their analysis. We repeat the description of their analysis from 
Branch et al. (2002) here. In their analysis, de Kool et al. (2002) used effective continuous 
spectra that consisted of a power-law continuum plus Fe II and Mg II Broad Emission Lines 
(BELs). The Mg II BEL was the sum of two Gaussians centered on the two components 
of the Mg II A2798A doublet (A2796A, A28031). Two different templates for the Fe II 
BELs were considered. The first consisted of a linear combination of five sets of Fe II 
BELs from theoretical model calculations (Verner et al. 1999), and the second was the 
observed Fe II BEL spectrum of the strong emission-line QSO 2226-3905 (Graham et al. 
1996). For the absorption features, de Kool et al. (2002) obtained a template distribution of 
line optical depth with respect to velocity in the Broad Absorption Line Region (BALR) from 
the observed absorption profile of Fe II A3004A, an apparently unblended line of moderate 
strength. Given the assumption that only absorption takes place in the BALR, the optical 
depth was obtained from r(v) = log F^, where F^ is the fractional residual flux in the 
absorption feature. The resulting optical depth distribution extended from about 1200 to 
2700 km s~ x and peaked near 2100 km s -1 . This optical depth distribution, scaled in 
amplitude, was used for all absorption lines. For each of the absorbing ions that were 
introduced Fe II, Mg II, Cr II, and Mn II the column density was a fitting parameter. The 
relative strengths of the lines of each ion turned out to be consistent with LTE. 

Three models that differed in their details (de Kool et al. 2002) were presented, with 
similar results. The column densities of Fe II, Cr II, and Mn II were well constrained. (The 
column density of Mg II could not be well constrained because the only Mg II absorption, 
due to A2798, is saturated.) The excitation temperature was found to be near 10,000 K. 
Two local covering factors, representing the fractions of the power-law source and the BELR 
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that are covered by the BALR as viewed by the observer, were introduced to reproduce 
the observed nonblack saturation, that is the fact that in the observed spectrum, even very 
strong absorption features do not go to zero flux. Both local covering factors were found to 
be 0.7 ±0.1. 

A detailed view of the spectral fit for one of the models was presented, and practically 
all of the observed absorptions were reasonably well accounted for. In order to interpret 
the results of FIRST J121442.3+280329. of their spectrum fits, de Kool et al. (2002) used 
the photoionization-equilibrium code Cloudy to compute a grid of constant-density slab 
models irradiated by a range of ionizing spectra. The ionization parameter U, the hydrogen 
density n and the hydrogen column density Nh were found to satisfy 2.0 < logU < 0.7, 
7.5 < logn < 9.5, and 21.4 < \ogN H < 22.2, respectively. From these values, the distance 
of the BALR from the center of the QSO was inferred to be between 1 and 30 pc. 

2.2.2. The PHOENIX approach 

PHOENIX calculations are very detailed, starting from first principles. We compare it with 
the rather simplified, template fitting approach of de Kool et al. (2002) (see also: Korista 
et al. 1992; Arav et al. 2001; de Kool et al. 2001) and discuss differences in the approaches. 
More details of PHOENIX are presented in an Appendix. 

In contrast to the approach by de Kool et al. (2002) we model the continuum, the 
emission lines, and the absorption lines simultaneously. The Fe II emission lines in our model 
are calculated self-consistently using our rather large model atom. The velocity width of the 
emission lines is determined by the outflow characteristics. In fact, the radiative transfer 
through the moving clouds is what gives the lines their width and therefore the calculations 
from the model are directly comparable with the observed spectrum. In addition we have 
run in excess of 100 models during the fitting process fitting the continuum, emission and 
absorption line features simultaneously, compared with the linear combination of the five 
calculations run with the sophisticated Cloudy models in the de Kool et al. (2002) approach. 
None of our models are based on observed emission line spectra and as such we do not have 
the problems associated with the correct placement and removal of the continuum. 

We have none of the fitting parameters associated with the emission line fits displayed in 
de Kool et al. (2002). We rely on a global fit to the overall spectrum from a model calculated 
from first principles, to determine how our model compares with observation. 

Our approach does not use a template fit to the absorption line spectrum r(v) = log F^. 
In our model the spectrum is calculated directly from self-consistently solving the radiative 
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transfer problem with scattering in an expanding atmosphere from first principles. The 
full spectral energy distribution, including continua, emission, absorption, and the effects of 
differential expansion are solved for simultaneously and then the results are compared to the 
observed spectra. Our optical depth profile as a function of velocity is determined based 
on the physical conditions of the gas and the radiation field at a given velocity, not on the 
observed spectrum of the object a priori. 

Our model grid is specified in terms of r st( j and the spatial extent is then determined 
by the implicit condition that T st ^(R )= 1.0. Thus, the inner and outer spatial boundaries 
are determined by the input r st d grid and vary not only from model to model but within a 
given model from iteration to iteration as we reach convergence. Given the density param- 
eterization, the pressure and temperature are determined by the NLTE equation of state 
from iteration to iteration. The models are very well converged. In contradistinction to the 
template fitting of method de Kool et al. (2002), since we do not have a spatially fixed inner 
boundary and the r st d grid is fixed, we do not have a priori, a fixed ionization parameter or 
column density. 



3. Determining the Best Fit Model 

In order to find the appropriate model parameters we began with the parameters of 
Branch et al. (2002) and restricted our initial calculations to pure LTE, which is much 
less computationally demanding and hence allows us to produce large grids which are then 
calculated in full NLTE. Given an LTE grid, we then chose a subset of "best fit" (by eye) 
models which we proceeded to calculate in full NLTE. Once we felt subjectively that we had 
the best model fit we turned to a more objective method of determining the best fit. Over 
100 models were run during this process. 

Determining a goodness of fit criterion for spectra is a non-trivial task. A pure x 2 
method is not ideal because a good fit of a synthetic spectrum to an observed spectrum 
should ideally fit on all scales. That is, the continuum shape (colors) should be correct as 
well as all the line features. Since our synthetic spectral models include detailed physics 
of the interaction of the lines and continua, we strive to fit both. Since we don't know the 
errors due to flux calibration, reddening, etc., a x 2 would not have the traditional meaning in 
terms of probability. Figure 3 illustrates how the direct comparison of models to observations 
clearly leads to the best fit spectrum. While the two models with T mo d e i= 4600 and T mo d c i= 
4700 look quite close, we used a x 2 -like test found the smallest distance between the model 
and the observed treating the normalization as a nuisance parameter. 
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4. Results 

4.1. FIRST J121442.3+280329 

We first turn our attention to the FeLoBAL quasar FIRST J121442.3+280329 (Becker 
et al. 2000a; White et al. 2003) UV and optical spectra. We adopted a redshift for FIRST J121442.3+280329 
of z — 0.692 as measured from the peak of the Mg II emission feature. We corrected for a 
Galactic extinction of E(B - V) = 0.023 (Schlegel et al. 1998) obtained from NED using 
the reddening curve of Cardelli et al. (1989) supplemented by the work of O'Donnell (1994). 
We smoothed all observed spectra using a near-Gaussian smoothing procedure with a width 
of 300 km s -1 . 

The model calculations and fits were done according to the method discussed in § 3. 
Figure 1 compares our best-fitting synthetic PHOENIX spectrum and the restframe UV spec- 
trum of FIRST J121442.3+280329. The best-fit parameters for FIRST J121442.3+280329 
were: T modcl = 4600 K, Rq = 1.4 x 10 17 cm, v e = 300 km s _1 , and v = 2100 km s _1 (for 
completeness these values are also shown in Table 1). 

With this radius and model temperature, the luminosity of the model is L = 6.3 x 
10 45 ergs s _1 . Overall, the synthetic model compares favorably with the restframe UV 
observation of FIRST J121442. 3+280329; in particular the Fe II lines, both from ground and 
excited states, fit very well. The Mg II A2798 feature appears to be too strong in emission and 
yet too weak in absorption. This could be a sign of asymmetry in the atmosphere, Branch 
et al. (2002) and de Kool et al. (2002) found a similar result in their calculations. Figure 2 
compares the FIRST J121442. 3+280329 rest-frame optical observations to the synthetic 
model. The synthetic model calculation is a reasonable fit to the optical observation. The 
features at 4600A are reproduced in the synthetic spectrum, yet are not as strong as observed 
and H I absorption appears to be too strong in the synthetic spectrum. However, this is 
dependent on our placement of the continuum and it should be noted that the observed 
optical spectrum is very noisy. The observations of Hall (2007) show that Balmer absorption 
features do exist in optical spectra of one FeLoBAL QSO which is very similar to the QSOs 
studied here. 

In order to reproduce the optical observation with our model we placed our lower bound- 
ary condition (that the specific intensity is given by the solution to the diffusion equation at 
the inner edge of our opaque core) at the somewhat low value of r st d = 10. Ideally we place 
that opaque core at r std ^ 100, but when we did that we found that the model spectrum 
longward of 2800 A no longer fit the observation. In particular the optical flux was very 
attenuated and had deep P-Cygni profiles unlike the observation. FeLoBALs, although very 
optically thick by AGN standards, are therefore thought to be optically thinner than, for 
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example, the atmosphere of a Type II supernova. Models in which the continuum optical 
depth is high {r std = 100) fail to reproduce the optical, and the UV spectra redward of Mg 
II of FeLoBALs. This could be an indication that our inner boundary condition would be 
better replaced by an AGN-like continuum, but that is beyond the scope of the present work. 

4-1.1. Importance of NLTE effects 

It is important to correctly model photoionization and recombination in the models. 
Figure 4 shows the importance of NLTE effects in correctly modeling spectra. The solid line 
has all the species that we included in these calculations in NLTE, whereas the dashed line 
has everything in NLTE except for Ca I-III. For species treated in LTE the Saha-Boltzmann 
equations are solved in order to calculate the atomic level populations instead of the full 
rate equations. An LTE treatment of calcium results in under-estimation of the ionization 
of calcium in the model which creates a large quantity of Ca II in the line forming region. 
Therefore the Ca II H&K features appear in the LTE spectrum. When NLTE is turned 
on for Ca I-III in the PHOENIX model, the level populations are controlled by the hotter 
radiation field (photoionized) from deeper layers and calcium is overionized compared to the 
local gas temperature. Therefore the Ca II H&K features, which are resonance transitions, 
disappear from the synthetic spectrum. The synthetic model which treats Ca I-III in NLTE 
closely reproduces the optical observation whereas the synthetic model without the NLTE 
Ca I-III clearly overestimates the strength of the H&K lines. 

4.2. ISO J005645. 1-273816 

The optical spectrum of ISO J005645. 1-2738 16 was obtained in September 2000 with 
the FORS1 instrument installed on the VLT UTl/Antu (Due et al. 2002). We adopted a 
redshift of z — 1.776 which was determined by the C IV, Fe II A2627 line, and Mg II emission 
lines. We corrected for an Galactic extinction of E(B — V) = 0.017 (Schlegel et al. 1998) 
using the standard reddening curve (Cardelli et al. 1989; O'Donnell 1994). We smoothed all 

Table 1. PHOENIX best-fitting model parameters. 



^0 v v e log(L bol ) 



1.4 x 10 17 cm 2100 km s" 1 300 km s" 1 45.8 ergs s" 1 
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observed spectra using a near-Gaussian smoothing procedure with a width of 300 km s l . 

The Mg II feature has a similar shape when compared with FIRST J121442.3+280329 
and ISO J005645. 1-2738 16 appears to have similar Fe II features when compared with 
FIRST J121442.3+280329. Because of the similarities between ISO J005645. 1-273816 and 
FIRST J121442.3+280329 we compare the synthetic model spectrum which fit FIRST J121442.3+280329 
with the rest frame UV spectrum of ISO J005645. 1-273816. This comparison is shown in 
Figure 5. Following the same fitting procedure outlined in § 3 and the same grid of models 
used for FIRST J121442. 3+280329 we found that the parameters which were a best fit for 
FIRST J121442.3+280329 also were a best fit for ISO J005645. 1-273816. The Fe II lines fit 
very well for this object and the Mg II emission feature was too strong while the absorption 
was too weak. This is very similar to the fit to the UV spectrum of FIRST J121442.3+280329. 

The similarities between the UV spectra of FIRST J121442.3+280329 and ISO J005645.1- 
273816 are interesting. FIRST J121442.3+280329 and ISO J005645.1-273816 may be a sub- 
type of FeLoBAL AGN with very similar characteristics. For emphasis, in Figure 6 we show 
the synthetic model spectrum with the combined FIRST J121442.3+280329 UV and optical 
spectra and the ISO J005645. 1-273816 UV spectrum. The high quality fit, over a wide wave- 
length range, is compelling. It is however unlikely that the PHOENIX models would compare 
favorably with the total composite spectrum of ISO J005645. 1-273816. ISO J005645. 1-273816 
was discovered in the infrared and the UV spectrum is diminished with respect to the in- 
frared (Due et al. 2002); however the IR emission could be enhanced by dust emission in a 
physically separate region. 

4.3. Physical Conditions 

The physical conditions in atmosphere are calculated from the solution of the radiative 
transfer equation, the equilibrium rate equations, and the equation of radiative equilibrium 
for the best fit input parameters given in Table 1. The PHOENIX model for ISO J005645.1- 
273816 and FIRST J121442. 3+280329 has the following physical dynamics: outflow mass 
550 Mq, kinetic energy 30 x 10 51 ergs, and a mass loss rate of M = 159 M Q yr _1 above the 
"photosphere" (r st d = 1)- In addition the PHOENIX model yields an outflow mass of 3000 M , 
kinetic energy of 100 x 10 51 ergs, and a mass loss rate of M = 466 M yr" 1 above T st <i = 10. 
The PHOENIX model has an equivalent hydrogen column density of 2 x 10 24 cm~ 2 for the 
region above the "photosphere" (r s td = 1) and a maximum equivalent hydrogen column 
density of 2 x 10 25 cm -2 for the entire atmosphere {r st d = 10). These values are displayed 
in Table 2. 
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5. Discussion 

The model fits are especially good considering the few free parameters in the PHOENIX 
model. The reader should keep in mind that we are not fitting each absorption line sepa- 
rately, but rather construct a global fit to the entire spectrum using the parameters given in 
Table 1. This includes the absorption, the emission, and the continuum, all fit with these few 
parameters. In fact, our model spectrum is created by the solution to the radiative transfer 
equation at every wavelength point across the spectrum simultaneously. 

In this section, we discuss some of the physical constraints from these models, how 
observed polarization in FeLoBALs fits in with this model, provide a few inferences about 
other FeLoBALs and then discuss some implications for quasar populations and evolution. 

5.1. Physical Constraints Inferred 

Although we have assumed a specific density profile, the radial extension of our model 
is very small (v max = 2800 km s -1 and v = 2100 km s -1 ); thus the model does not strongly 
probe the density structure of the ejecta, except to infer that the density profile is rather 
flat. Branch et al. (2002) assumed a power-law density profile with p oc (^)~ n , with n = 
2, whereas our effective power-law index at the pseudo-photosphere r st d = 1 is vo/v e = 
7, nominally significantly steeper, but since the radial extension of our model is so small 
the discrepancy should not be considered important. What is perhaps more interesting is 
that our value of v is quite a bit higher than that of Branch et al. (2002) who used a 
photospheric velocity of 1000 km s" 1 . Interestingly, they had to impose a minimum velocity 
of 1800 km s _1 on Fe II and Cr II and they chose the same value of v max = 2800 km s _1 . 
The lower photospheric velocity is likely to be a consequence of the Schuster-Schwarzschild 
approximation (Mihalas 1978) of SYNOW. That is, in SYNOW all emission of photons occurs 
at the photosphere and the atmosphere consists merely of a "reversing-layer" where there is 
no creation or destruction of photons, only resonant scattering. In contrast, PHOENIX allows 



Table 2. PHOENIX Physical Conditions. 



r std 


Outflow Mass 


Kinetic Energy 


M 


Column Density 


1 


547 M 


30 x 10 51 ergs 


159 M yr' 1 


2 x 10 24 cm" 2 


10 


3000 M 


100 x 10 51 ergs 


466 M Q yr" 1 


2 x 10 25 cm" 2 
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line formation throughout the atmosphere, and lines typically form in the "line-forming 
region" (3.0 < r stc j < 0.1) depending on their strength. The much lower "photospheric" 
temperature that we find is robust since SYNOW does no continuum transfer and thus 
the "temperature" in a SYNOW model is not physically meaningful; it just is a way of 
parameterizing the underlying continuum. However, PHOENIX solves the full NLTE radiative 
transfer problem and particularly our result that the Ca II H+K feature is only reproduced 
in NLTE indicates that we have the right physical conditions in our model. 

The synthetic model has 547M and kinetic energy 30 x 10 51 ergs above the "photo- 
sphere" (r s td = 1)- With i?o — 1-4 x 10 17 cm and maximum velocity v max = 2800 km s" 1 
we estimate a crossing time of t ~ Ro/ v max = 15.5 yr. We estimate the mass loss rate 
using M ~ M/t = 35 M Q yr" 1 which is 1/5 the mass loss rate in our PHOENIX model 
of M — 159 M yr" 1 . The kinetic energy luminosity is estimated to be E k ~ E k /t = 
4.3 x 10 43 erg s" 1 which is two orders of magnitude lower than the bolometric luminosity of 
the model L^i = 6.3 x 10 45 erg s _1 ; thus the flow could be luminosity driven. Even more 
interesting is that the velocity we find is very similar to characteristic velocities of hot stellar 
winds which are thought to be driven by line absorption in the atmosphere (Walborn et al. 
1995). 

Branch et al. (2002) found a luminosity of 6 x 10 46 erg s _1 from the photometry of 
FIRST J121442. 3+280329 and using the quasar composite spectrum of Mathews & Ferland 
(1987) to perform the K-correction and the relation L^i = 9AL^ at 5100A (Kaspi et al. 2000). 
Using the quasar composite spectrum of Francis et al. (1991) we find Lboi = 4.4 x 10 46 erg s" 1 
for FIRST J121442.3+280329 and L bo i = 1.7 x 10 45 erg s" 1 for ISO J005645. 1-273816. Since 
the spectra are so similar, but the total luminosities inferred differ by a factor of 10, sug- 
gesting that the K-correction is important, as it is dependent on the shape of the spectral 
energy distribution. 

Our models have a maximum column density of 2 x 10 25 cm" 2 for the entire atmosphere 
and a column density of 2 x 10 24 cm" 2 for the region above the "photosphere" r st d = 1- 
As discussed in § 4.1, the maximum column density is constrained by our fit to the optical 
spectrum. A higher column density would yield higher continuum and line optical depths 
and the continuum would be much weaker than observed. Thus, the requirement that a 
photosphere be present places a limit on the minimum column density of 2 x 10 24 cm -2 , and 
the optical spectrum places a limit on the maximum column density of 2 x 10 25 cm" 2 of the 
wind. These constraints make these models Thompson thick, and is sufficient to make these 
objects appear to be X-ray faint as observed (Green et al. 2001; Gallagher et al. 2006). 

We can now calculate a luminosity distance using a variant of the Spectral-fitting Ex- 
panding Atmosphere Method (SEAM) used to derive distances to supernovae (Baron et al. 
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1995; Lentz et al. 2001a; Baron et al. 1999; Mitchell et al. 2002; Baron et al. 2004). SEAM 
is a sophisticated variant of the classical Baade-Wesselink method (Baade 1926). Our ap- 
proach here differs from the conventional SEAM method in the crucial respect that in the 
supernova case, we know both that homology is an excellent approximation shortly after 
the explosion and we know that there was an explosive ejection at a time to (which SEAM 
determines). Here, homology has been taken to be an expedient Ansatz and Rq has been set 
semi-empirically. There is no reason a priori to expect that FeLoBALs are the result of a 
single ejection event and thus the radius in our models is much more poorly known than in 
the case of supernovae. However the radial extent does determine the overall density, and if 
we were able to identify features which are good density indicators we may be able to place 
the results of this method on firmer ground. Of course, we are sensitive to systematic errors 
in the overall flux calibration and in the total reddening. Errors due to reddening tend to 
cancel out, since the more we must deredden the observed spectrum, the hotter we must 
make the synthetic spectrum, which compensates for the dimmer observed spectrum. Thus, 
below we will test the value used for R by using the synthetic spectra to calculate synthetic 
photometry and the predicted bolometric magnitudes in a number of bands. Comparing 
these magnitudes to the observed photometry, we obtain a distance. This result is sensitive 
to Ro and thus if we find reasonable values for the distance we can have confidence in our 
choice of Ro, and hence that we have the right density. The outflow rate we find is then also 
reasonable. Thus, this provides an important check on our parameters. 

The SEAM method uses observed photometry and the synthetic spectrum to calculate 
synthetic photometry as well as K-corrections. From this we find a distance modulus ji which 
is simply related to the luminosity distance. Using m B = 17.06 for FIRST J121442.3+280329 
we find \x = 42.56 or di = 3.25 Gpc, which compares favorably with the luminosity distance 
inferred for our adopted cosmology (H = 70, Qm — 0.3, Vt^ — 0.7) of = 4.2 Gpc. Using 
m B = 22.74 for ISO J005645.1-273816 we find /i = 46.59 or d L = 20.8 Gpc, which is a bit high 
compared with the luminosity distance inferred for our adopted cosmology d^ = 13.4 Gpc. 
As we noted above, it seems likely that the K-corrections are important and thus we obtain 
a distance for the relatively nearby FIRST J121442.3+280329 which is good to about 30%, 
but for the more distant ISO J005645. 1-273816 the distance is off by 55%. Nevertheless, that 
fact that the distances agree to better than a factor of two indicates that our model predicts 
roughly the right size as well as the right SED for both objects. 

Since our PHOENIX synthetic spectra are quite a good fit, we may also determine the 
bolometric luminosity using our synthetic spectra to perform synthetic photometry and 
calculate K-corrections. With this method we find that L^i = 1.2 x 10 46 erg s _1 for 
FIRST J121442.3+280329 and L bol = 2.9 x 10 45 erg s" 1 for ISO J005645.1-273816. 
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We obtained Lboi in two different ways, both of which were consistent with each other. 
However the two objects have nearly identical SEDs in the regions that we can observe 
(Fig. 6). Why does the synthetic L^i differ by a factor of 10? The values for R which are 
the same for the synthetic models must actually differ slightly between the two objects, or, 
as noted above, the K-correction is important. 

We use our model luminosity to estimate a few quantities. If we assume that our model 
luminosity is the Eddington luminosity then the black hole mass and the accretion rate are 
Me = 5 x 10 7 M and Me = 1.1 M yr _1 . These are roughly consistent with usual estimates 
for quasars, even if we scale the luminosity up by a factor of 3 that we infer from photometry. 

5.2. Polarization 

Lamy & Hutsemekers (2004; 2000) found that the polarization increases in the ab- 
sorption troughs do not rule out the model described in this paper. They constructed a 
"two-component" wind model. In their model, the broad absorption occurs in a dense equa- 
torial wind emerging from the accretion disk, while scattering and polarization mainly take 
place in a polar region. Our model is consistent with the two component model (Lamy & 
Hutsemekers 2004) in which the observer looks through a equilateral wind and sees a polar 
component which is dominated by electron scattering. A spherically symmetric model is a 
required computational constraint inherent in PHOENIX; thus we assume 100% global cover- 
ing. The polarization results indicate that some asymmetry must be present; nevertheless, 
the presence of P-Cygni profiles where the absorption troughs go almost to zero flux indicate 
that the covering fraction is high. 

Our model will not change much if we relax the 100% global covering fraction and have 
the same electron scattering polar component described in Lamy & Hutsemekers (2004). In 
fact, reducing the covering factor would most likely provide a better fit, as it would reduce 
the emission feature in Mg II that was shown in §4. 1 to be slightly too big. 

5.3. Other FeLoBAL QSOs 

FeLoBALs are observed to have a wide range of spectra (Hall et al. 2002). We fit only 
the spectra of two objects in this paper; in this section, we comment briefly on whether or 
not the PHOENIX model is likely to be able to explain the spectra of other objects. 

Hall (2007) report blueshifted broad absorption lines troughs in Balmer lines in the 
quasar SDSS J125942.80+121312.6. Our PHOENIX model predicts Balmer absorption in at 
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least H/3 and H7 as seen in Fig. 2. Thus our model may very well be able to explain the 
spectra of this object too. 

The value of the global covering fraction is very likely to be responsible for the differences 
between the spectra that we have modelled here, and those of other FeLoBAL QSOs (Hazard 
et al. 1987; Cowie et al. 1994; Becker et al. 1997; Hall et al. 2002). Specifically many 
FeLoBALs appear to have only absorption features, and the emission features are week. 
The high global covering fraction inherent in the PHOENIX model predicts rather prominent 
emission features and P-Cygni profiles. Reducing the global covering while viewing the 
outflow from the equatorial region would reduce the emission line strengths while retaining 
the strong absorption features. 

So-called overlapping trough QSOs (Hall et al. 2002) may also be able to be explained 
by this paradigm. We investigated them with our models and yet where unable to come up 
with a satisfactory solution, because the emission in our 100% global covering method is too 
high. Therefore the covering fraction in these type of quasars might be not quite as large 
as for the objects modeled here. In addition they appear to have higher velocity winds and 
possibly higher optical depths. 

If the covering fraction is not 100%, then from some lines of sight we may see the nucleus 
directly. The strong Fe II absorption observed when the outflow is in the line of sight may 
be seen as as emission in this case, and ultrastrong Fe II emission may be seen. The idea of 
linking Fe II emission with BAL outflows not in the line of sight has been suggested for the 
prototypical Narrow-line Seyfert 1 (NLS1) I Zw 1 (Boroson & Meyers 1992). 

5.4. Implications for Quasar Populations and Evolution 

Broad absorption lines are found in the spectra of 10-20% of the quasars in the Sloan 
Digital Sky Survey (e.g. Trump et al. 2006). The simplest interpretation of this fact is that 
all quasars have an outflow that occults 10-20% of quasar lines of sight (the orientation 
model). Alternatively, the outflow may cover a larger percentage of sight lines, and the 
broad absorption lines may only be present during some small fraction of the quasar lifetime 
during which it is blowing gas out of the nucleus (e.g., Voit et al. 1993; Becker et al. 2000b; 
Gregg et al. 2006). As discussed in § 5.2, our model requires a large global covering fraction, 
although it does not have to be 100%. Here we briefly review the support for both scenarios, 
and discuss how our results fit into models for quasars in general. 

Weymann et al. (1991) compare the emission-line and continuum properties of spectra 
from 42 BALQSOs with those of 29 normal QSOs. They find that the emission-line and 
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continuum properties are very similar between the BALQSOs and the non-BALQSOs. This 
result forms a key support for the orientation model. Another important piece of evidence 
that supports the idea that BALQSOs and non-BALQSOs differ only in their orientation 
is the fact that X-ray spectra from some BALQSOs are highly absorbed, but apparently 
intrinsically identical to those of non-BALQSOs (Gallagher et al. 2002). 

However, the Weymann et al. (1991) sample is very small, and the number of BALQSOs 
for which X-ray spectra of sufficient quality to obtain absorption column information is also 
small. Thus, while the orientation model is widely accepted, and may be applicable to 
many BALQSOs, the evidence for the evolutionary model, at least for some subclasses of 
BALQSOs, is growing. The best candidate for objects characterizing the evolutionary model 
and not fitting into the orientation model are the low-ionization BALQSOs. Early on, these 
were noted to have optical and UV emission-line and continuum properties different than 
ordinary quasars (Weymann et al. 1991; Boroson & Meyers 1992); for example, they tend 
to have strong Fe II emission and weak O III]. They also tend to be more reddened than 
non-BALQSOs and HiBALs (Reichard et al. 2003), and they are uniformly more X-ray weak 
than HiBALs (Gallagher et al. 2006). In addition, they are very rare; Trump et al. (2006) 
find only about 1% of quasars are LoBALs. Thus, the usual argument used in favor of 
the orientation model, the similarity of spectra from BALQSOs and non-BALQSOs, doesn't 
work as well for LoBALs. 

Further evidence for an evolutionary role of BALQSOs comes from their radio properties. 
If BALQSOs were observed predominately edge-on, as are radio galaxies, one would expect 
to see a steep radio spectrum dominated by the synchrotron emission in the radio lobes. 
However, BALQSOs show both steep and flat radio spectra (e.g., Becker et al. 2000b). 
Furthermore, Fe II BALQSOs are extremely rare, and there is evidence among the small 
sample of an anti-correlation between the radio-loudness and the strength of the BAL features 
that led Gregg et al. (2006) to propose that quasars in this state are emerging from cocoons 
of gas that produces the BALs and which suppresses the development of radio jets and lobes. 
In addition, Brotherton et al. (2006) present spectropolometric results showing polarization 
parallel to the radio axis, implying a small angle of inclination; they present an extensive 
review of the implications of the radio properties on BALQSO models. 

Numerous models of quasar evolution admit a time period early in the life of a quasar 
when it is heavily shrouded by dust and gas. Before we see the bare quasar, we expect to see 
it in a heavily absorbed stage. Outflows may result from the turning on of the QSO (Hazard 
et al. 1984). LoBALs may be young quasars that are casting off their cocoons of dust and 
gas (Voit et al. 1993) and may be related to ultra-luminous infrared galaxies (Egami 1999). 
Becker et al. (1997) suggest that FeLoBALs may be the missing link between galaxies and 
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quasars. Recent simulations of mergers and quasar evolution show that during much of the 
lifetime of the quasar, it is heavily shrouded by gas with large column densities (e.g., Hopkins 
et al. 2005). Observational support for this view comes from the discovery of a large number 
(four out of a sample of eight) LoBALs in high redshift quasars (Maiolino et al. 2004). 

Thus, observational evidence and scenarios for quasar evolution imply that FeLoBALs 
may be candidates for support of the evolutionary scenario. If the quasar observable lifetime 
is assumed to be 10 7 years, the fact that FeLoBALs only make up 0.33% of quasars in the 
SDSS (Trump et al. 2006) implies that the time scale for this stage is approximately 3 x 10 4 
years. That should be a lower limit, however, since these objects may be missed in the SDSS 
(Trump et al. 2006). This relatively short time period, compared with the total lifetime of 
a quasar, compares favorably with our crossing time of 15 years. 

It is clear that a high global covering fraction is compatible with the evolutionary sce- 
nario. But what about the other parameters obtained from the model fits?. One possible 
problem is that the radius of the model shows that the outflow occurs well within the central 
engine, at approximately 1.4 x 10 17 cm. The quasar feedback scenario proposed by Fabian 
(1999) infers that the bulk of the gas expelled from the AGN is accelerated beyond the Bondi 
radius, which for a 5 x 10 7 M Q black hole is 4 x 10 21 cm. On the other hand, the amount 
of gas expelled during this process could be comparable to the amount of gas contained by 
the entire galaxy, far more than the rather modest 547M Q that we infer. It is quite possible 
that the accretion/blow out process is messy and chaotic; there may be a phase in which the 
central engine "burps" , and this results in the features that we see. 

6. Conclusions 

Using the spectral synthesis code PHOENIX we compute synthetic spectra that provide 
very good fits to the observed restframe UV and optical spectra of two FeLoBALs. While 
our models are limited to exact spherical symmetry, they provide excellent fits. In order to 
reconcile our results with the polarization data on these objects, we require some asymmetry, 
but still a high global covering fraction, which would only modestly affect the flux spectrum. 

We are able to determine a luminosity distance estimate which is direct and is accurate 
to around 50%. The question arises: could these objects be used as distance indicators at 
high z, even if only as a sanity check on the really high-z Hubble diagram from GRBs? 

Our results lend support to the inference that FeLoBALs are an evolutionary stage of 
the QSO as opposed to a pure orientation effect. Our model with a smaller covering factor 
may be able to explain other BAL QSO such as overlapping trough QSOs. In addition our 
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model column densities, which are Compton thick, match those that are expected from X-ray 
non-detections of these objects. 

For future work we plan to explore metallicity effects on the model spectrum. This 
is not trivial, as it is not just a matter of changing the metallicity and comparing the 
model. Completely different sets of dynamical and luminosity parameters may be required 
to achieve the best fitting model spectrum. We plan to continue our analysis with the Hall 
(2007) spectrum, which shows Hf3 absorption. We have compared our observed spectrum 
with the one from that paper and think that it is an excellent candidate for this type of 
modeling. 
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Fig. 1. — The PHOENIX spectrum (solid line) compared with the restframe, dereddened, and 
smoothed, observed FIRST J121442.3+280329 (dotted line) UV spectrum. 




Fig. 2. — The PHOENIX model (solid line) compared with the restframe, dereddened, and 
smoothed, observed optical spectrum of FIRST J121442.3+280329 (dotted line). 
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Fig. 3. — Four models (black lines) are shown, normalized to the observed spectrum of 
FIRST J1214+2803 (grey lines) at 3600A. The top and bottom models are seen to provide 
a distinctly poorer fit than the middle two models, especially in the 2300-2900A band. The 
T = 4600 K and T = 4700 K provide a much better fit, and although they appear very 
similar, the T = 4600 K was selected as the best fit by minimizing a x 2 -like figure of merit. 
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Fig. 4.— The optical spectrum of FIRST J121442.3+280329 (dotted) deredshifted, dered- 
dened, and smoothed, versus two synthetic PHOENIX models. The dashed line represents a 
PHOENIX model with the ions discussed in §2 in NLTE except for the Calcium I-III ions, which 
are in LTE. The solid line shows a model with all of the ions in NLTE that are discussed in 
§2. The fact that the synthetic spectrum matches the observed spectrum significantly better 
indicates the importance of NLTE in the calculations. 
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Fig. 5. — The PHOENIX model (solid line) vs the restframe, dereddened, and smoothed, UV 
spectrum of ISO J005645.1-273816 (dotted line). 
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Fig. 6.— The combined spectra of FIRST J121442.3+280329 and ISO J005645. 1-273816. 
The wavelength range spans 1500-5500A. the PHOENIX spectrum is the solid line, while the 
FIRST J121442.3+280329 UV through optical spectrum and the ISO J005645.1-273816 UV 
spectrum are dotted. All spectra are restframe, deredshifted, and smoothed. 
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A. The PHOENIX Code 

PHOENIX is a mature code. Development of PHOENIX began in 1990 and continues 
today. Discussions of the computational details appear in a number of publications. We 
include a short description of some of the computational details for those who are interested 
in understanding in more detail how the code works. This section also includes relevant 
references. Finally, we present a flowchart of a PHOENIX computation. 

A.l. An introduction to PHOENIX 

PHOENIX is able to model astrophysical plasmas under a variety of conditions, in- 
cluding differential expansion at relativistic velocities (Hauschildt & Baron 2006, 2004a; 
Baron & Hauschildt 2004; Hauschildt & Baron 1999; Hauschildt et al. 1997a,b; Allard et al. 
1997). PHOENIX includes very detailed model atoms constructed from the work of Kurucz 
(Hauschildt & Baron 1995; Baron et al. 1997; Short et al. 1999) for a number of almost all 
the important species (H, He, CNO, Si, Fe, Co, etc.). In addition, the CHIANTI or APED 
databases may be chosen for model atoms at runtime. The code is optimized and parallelized 
to run on all available supercomputers. PHOENIX has a long history of modeling astrophysical 
objects including extra-solar giant planets (EGPs), Brown dwarfs (Allard et al. 1997, 2001), 
novae (Petz et al. 2005), as well as all types of supernovae (Baron et al. 1995, 1999; Nugent 
et al. 1997; Baron et al. 2000; Lentz et al. 2001a,b,c; Mitchell et al. 2001; Baron et al. 2006). 
In version 14, we solve the fully relativistic radiative transport equation for a variety of spa- 
tial boundary conditions in both spherical and plane-parallel geometries for both continuum 
and line radiation simultaneously and self-consistently using an operator splitting technique. 
We also use an operator splitting technique to solve the full multi-level NLTE transfer and 
rate equations for a large number of atomic species (with a total of more than 10,000 energy 
levels and more than 100,000 primary NLTE lines), including non-thermal processes. MPI 
and OpenMP directives are used, so the code runs on both distributed and shared memory 
architectures (Baron et al. 2003; Hauschildt et al. 2001; Baron & Hauschildt 1998; Hauschildt 
et al. 1997a). PHOENIX accurately solves the fully relativistic radiation transport equation 
along with the non-LTE rate equations (currently for ~ 150 ions) while ensuring radiative 
equilibrium (energy conservation). Typically each atom has several ionic species in NLTE 
and is represented by dozens to hundreds of levels for the Fe-group species. PHOENIX is 
currently around 700,000 lines of code which relies on 0.6 GB of atomic data and 12 GB of 
molecular data. 

In the present paper, the multilevel, non-LTE rate equations are solved self-consistently 
for H I, He I— II, Mg I — III, Ca I-III, and Fe I III using an accelerated lambda iteration 
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(ALI) method (Rybicki & Hummer 1991; Hauschildt 1992; Hauschildt & Baron 1999, 2004b). 
Simultaneously we solve for the special relativistic condition of radiative equilibrium (Nugent 
et al. 1997) using a modified Unsold-Lucy temperature correction scheme. Relativistic effects, 
in particular the effects of advection and aberration, are important in the high velocity flows 
observed in these quasars. 

The generalized non-LTE equation of state (EOS) is solved for 40 elements and up to 26 
ionization stages per element for a total of hundreds of species. For the conditions present 
in the models, molecules are unimportant, and we neglect them in order reap substantial 
savings in CPU time. Negative ions are always included. The numerical solution of the EOS 
is based on Brent's method for the solution of nonlinear equations which is very robust and 
fast. 

In addition to the non-LTE lines, the models include, self-consistently, line blanketing 
of the most important (~ 10 6 ) lines selected from the latest atomic and ionic line list of 
Kurucz. The entire list contains close to 42 million lines but not all of them are important 
for the case at hand. Therefore, before every temperature iteration, a smaller list is formed 
from the original list. A set of optical depths in the line- forming region of the gas is chosen, 
then using the density and temperature for these depths, the absorption coefficient in the 
line center, Ki, is calculated for every line and compared to the corresponding continuum 
(LTE+NLTE) absorption coefficient, k c . A line is transferred to the "small list" if the ratio 
k,i/k c is larger than a pre-specified value (in these calculations 5 x 10 -6 , selecting over half a 
million lines). In the subsequent radiative transfer calculations all lines selected in this way 
are taken into account as individual lines and all others from the large line list are neglected. 
This selection procedure is repeated at each iteration where the pressure or temperature 
changes by a prescribed amount in order to always include the most important lines. We 
treat line scattering in these LTE lines by setting the albedo for single scattering, a = 0.95. 



A. 2. A PHOENIX Flowchart 

This section gives the reader a brief explanation of how the PHOENIX code computes 
a model spectrum. For a more complete understanding the authors recommend the reader 
peruse Hauschildt & Baron (1999). Our iteration scheme for the solution of the multi- level 
non-LTE problem can be summarized as follows: (1) for given population levels [n*] and elec- 
tron densities [n e ], solve the radiative transfer equation at each wavelength point and update 
the radiative rates and the approximate rate operator, (2) solve the linear system for the 
atomic level populations for a given electron density, (3) compute new electron densities by 
the generalized partition function method, (4) calculate the temperature corrections needed 
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to bring the current iteration into radiative equilibrium, (5) repeat until a fixed number of 
iterations is reached (and check that all quantities have converged). It is crucial to account 
for coherent scattering processes during the solution of the wavelength dependent radiative 
transfer equation, it explicitly removes a global coupling from the iterations. 

As the first step in our outermost iteration loop (the model iteration) we use the current 
best guess of [T, n,] as function of radius to solve the hydrostatic or hydrodynamic equations 
to calculate an improved run of P gas with radius. Simultaneously, the population numbers 
are updated to account for changes in P gas . The next major step is the computation of the 
radiation field for each wavelength point (the wavelength loop), which has the prerequisite 
of a spectral line selection procedure for LTE background lines. Immediately after the 
radiation field at any given wavelength is known, the radiative rates and the rate operators 
are updated so that their calculation is finished after the last wavelength point. In the next 
steps, the population numbers are updated by solving the rate equations for each NLTE 
species and new electron densities are computed, this gives improved estimates for [n«] . The 
last part of the model iteration is the temperature correction scheme outlined above (using 
opacity averages etc. that were computed in the wavelength loop) which delivers an improved 
temperature structure. If the errors in the constraint equations are larger than a prescribed 
accuracy, the improved [T, rij] are used in another model iteration. Using this scheme, about 
10-20 model iterations are typically required to reach convergence to better than about 1% 
relative errors, depending on the quality of the initial guess of the independent variables and 
the complexity of the model. 
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Flowchart for a global iteration of PHOENIX. 



